library(rdrobust)
library(rdd)
library(tidyverse)
library(estimatr)
library(modelsummary)
library(kableExtra)

load("RD_parties_and_violence.RData")


###################################################################
###################################################################
#       Why Programmatic Parties Reduce Criminal Violence:
#                 Theory and Evidence from Brazil
#
#                      Supplementary material
#
#          Camilo Nieto-Matiz                Natán Skigin
#       camilo.nieto-matiz@utsa.edu          nskigin@nd.edu
###################################################################
###################################################################


# code to export RD models with msummary()
tidy.rdrobust <- function(model, ...) {
  ret <- data.frame(
    term = row.names(model$coef),
    estimate = model$coef[, 1],
    std.error = model$se[, 1],
    p.value = model$pv[, 1]
  )
  row.names(ret) <- NULL
  ret
}

glance.rdrobust <- function(model, ...) {
  ret <- data.frame(
    Kernel = model$kernel,
    Bandwidth = model$bwselect
  )
  ret
}

#================================================
#                 Figure B.1
#              Covariate balance
#================================================

# local linear specification
summary(bal.test1<-rdrobust(y = bra_all$frac_landless, x = bra_all$prog, 
                            c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                            all = T, covs = NULL)) 
summary(bal.test2<-rdrobust(y = bra_all$banks, x = bra_all$prog, 
                            c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                            all = T, covs = NULL)) 
summary(bal.test3<-rdrobust(y = bra_all$unemp, x = bra_all$prog, 
                            c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                            all = T, covs = NULL)) 
summary(bal.test4<-rdrobust(y = bra_all$land_gini, x = bra_all$prog, 
                            c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                            all = T, covs = NULL))
summary(bal.test5<-rdrobust(y = bra_all$gini, x = bra_all$prog, 
                            c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                            all = T, covs = NULL))
summary(bal.test6<-rdrobust(y = bra_all$rain_yearly, x = bra_all$prog, 
                            c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                            all = T, covs = NULL)) 
summary(bal.test7<-rdrobust(y = bra_all$hdi_education, x = bra_all$prog, 
                            c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                            all = T, covs = NULL)) 
summary(bal.test8<-rdrobust(y = bra_all$log_security_budget, x = bra_all$prog, 
                            c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                            all = T, covs = NULL)) 
summary(bal.test9<-rdrobust(y = bra_all$log_social_budget, x = bra_all$prog, 
                            c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                            all = T, covs = NULL)) 
summary(bal.test10<-rdrobust(y = bra_all$ag_income, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test11<-rdrobust(y = bra_all$share_pib_nonag, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test12<-rdrobust(y = bra_all$frac_ext_pov, x = bra_all$prog, 
                             c = 0, p=1, vce = "nn", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test13<-rdrobust(y = bra_all$distcap, x = bra_all$prog, 
                             c = 0, p=1, vce = "nn", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test14<-rdrobust(y = bra_all$invasions_count, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test15<-rdrobust(y = bra_all$reforms_count, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test16<-rdrobust(y = bra_all$numfarmsratio_1_100ha, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test17<-rdrobust(y = bra_all$ruralpc, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL))
summary(bal.test18<-rdrobust(y = bra_all$muniguard_dum_min, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL))
summary(bal.test19<-rdrobust(y = bra_all$lnagprodareacult, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test20<-rdrobust(y = bra_all$cattle_depend, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test23<-rdrobust(y = bra_all$coffee_depend, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test24<-rdrobust(y = bra_all$leftgovernor, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test25<-rdrobust(y = bra_all$rightgovernor, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test26<-rdrobust(y = bra_all$psdb.vs.1994, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 
summary(bal.test27<-rdrobust(y = bra_all$pt.vs.1994, x = bra_all$prog, 
                             c = 0, p=1, vce = "hc0", bwselect = "mserd", kernel="tri",
                             all = T, covs = NULL)) 




#ols - full sample
bra_all$tr<-ifelse(bra_all$prog>0,1,0)
summary(bal.test1.f<-lm_robust(frac_landless~prog*tr, data=bra_all)) 
summary(bal.test2.f<-lm_robust(banks~prog*tr, data=bra_all)) 
summary(bal.test3.f<-lm_robust(unemp~prog*tr, data=bra_all)) 
summary(bal.test4.f<-lm_robust(land_gini~prog*tr, data=bra_all))
summary(bal.test5.f<-lm_robust(gini~prog*tr, data=bra_all))
summary(bal.test6.f<-lm_robust(rain_yearly~prog*tr, data=bra_all)) 
summary(bal.test7.f<-lm_robust(hdi_education~prog*tr, data=bra_all)) 
summary(bal.test8.f<-lm_robust(log_security_budget~prog*tr, data=bra_all)) 
summary(bal.test9.f<-lm_robust(log_social_budget~prog*tr, data=bra_all)) 
summary(bal.test10.f<-lm_robust(ag_income~prog*tr, data=bra_all)) 
summary(bal.test11.f<-lm_robust(share_pib_nonag~prog*tr, data=bra_all)) 
summary(bal.test12.f<-lm_robust(frac_ext_pov~prog*tr, data=bra_all)) 
summary(bal.test13.f<-lm_robust(distcap~prog*tr, data=bra_all)) 
summary(bal.test14.f<-lm_robust(invasions_count~prog*tr, data=bra_all)) 
summary(bal.test15.f<-lm_robust(reforms_count~prog*tr, data=bra_all)) 
summary(bal.test16.f<-lm_robust(numfarmsratio_1_100ha~prog*tr, data=bra_all)) 
summary(bal.test17.f<-lm_robust(ruralpc~prog*tr, data=bra_all)) 
summary(bal.test18.f<-lm_robust(muniguard_dum_min~prog*tr, data=bra_all)) 
summary(bal.test19.f<-lm_robust(lnagprodareacult~prog*tr, data=bra_all)) 
summary(bal.test20.f<-lm_robust(cattle_depend~prog*tr, data=bra_all)) 
summary(bal.test23.f<-lm_robust(coffee_depend~prog*tr, data=bra_all)) 
summary(bal.test24.f<-lm_robust(leftgovernor~prog*tr, data=bra_all)) 
summary(bal.test25.f<-lm_robust(rightgovernor~prog*tr, data=bra_all)) 
summary(bal.test26.f<-lm_robust(psdb.vs.1994~prog*tr, data=bra_all)) 
summary(bal.test27.f<-lm_robust(pt.vs.1994~prog*tr, data=bra_all)) 






#frame
balance <- data.frame(Var =   c(bal.test1.f$outcome,bal.test2.f$outcome,bal.test3.f$outcome,bal.test4.f$outcome,bal.test5.f$outcome,
                                bal.test6.f$outcome,bal.test7.f$outcome,bal.test8.f$outcome,bal.test9.f$outcome,bal.test10.f$outcome,
                                bal.test11.f$outcome,bal.test12.f$outcome,  bal.test13.f$outcome, bal.test14.f$outcome,
                                bal.test15.f$outcome, bal.test16.f$outcome, bal.test17.f$outcome, bal.test18.f$outcome,
                                bal.test19.f$outcome,  bal.test20.f$outcome, 
                                bal.test23.f$outcome, bal.test24.f$outcome, bal.test25.f$outcome, bal.test26.f$outcome,
                                bal.test27.f$outcome),
                      Estimate =   c(bal.test1$coef[3],bal.test2$coef[3],bal.test3$coef[3],bal.test4$coef[3],bal.test5$coef[3],
                                     bal.test6$coef[3],bal.test7$coef[3],bal.test8$coef[3],bal.test9$coef[3],bal.test10$coef[3],
                                     bal.test11$coef[3],bal.test12$coef[3],bal.test13$coef[3],bal.test14$coef[3],bal.test15$coef[3],
                                     bal.test16$coef[3],bal.test17$coef[3],bal.test18$coef[3],bal.test19$coef[3],bal.test20$coef[3],
                                     bal.test23$coef[3],bal.test24$coef[3],bal.test25$coef[3],
                                     bal.test26$coef[3],bal.test27$coef[3]),
                      pvalue =     c(bal.test1$pv[3],bal.test2$pv[3],bal.test3$pv[3],bal.test4$pv[3],bal.test5$pv[3],
                                     bal.test6$pv[3],bal.test7$pv[3],bal.test8$pv[3],bal.test9$pv[3],bal.test10$pv[3],
                                     bal.test11$pv[3],bal.test12$pv[3],bal.test13$pv[3],bal.test14$pv[3],bal.test15$pv[3],
                                     bal.test16$pv[3],bal.test17$pv[3],bal.test18$pv[3],bal.test19$pv[3],bal.test20$pv[3],
                                     bal.test23$pv[3],bal.test24$pv[3],bal.test25$pv[3],
                                     bal.test26$pv[3],bal.test27$pv[3]),
                      SE =         c(bal.test1$se[3],bal.test2$se[3],bal.test3$se[3],bal.test4$se[3],bal.test5$se[3],
                                     bal.test6$se[3],bal.test7$se[3],bal.test8$se[3],bal.test9$se[3],bal.test10$se[3],
                                     bal.test11$se[3],bal.test12$se[3],bal.test13$se[3],bal.test14$se[3],bal.test15$se[3],
                                     bal.test16$se[3],bal.test17$se[3],bal.test18$se[3],bal.test19$se[3],bal.test20$se[3],
                                     bal.test23$se[3],bal.test24$se[3],bal.test25$se[3],
                                     bal.test26$se[3],bal.test27$se[3]))
balance$Coefficient <- ifelse(balance$Estimate>0, "Positive", "Negative")
balance$Var <- factor(balance$Var, ordered(balance$Var))
balance$Specification <- "Local"

balance1 <- data.frame(Var =   c(bal.test1.f$outcome,bal.test2.f$outcome,bal.test3.f$outcome,bal.test4.f$outcome,bal.test5.f$outcome,
                                 bal.test6.f$outcome,bal.test7.f$outcome,bal.test8.f$outcome,bal.test9.f$outcome,bal.test10.f$outcome,
                                 bal.test11.f$outcome,bal.test12.f$outcome,  bal.test13.f$outcome, bal.test14.f$outcome,
                                 bal.test15.f$outcome, bal.test16.f$outcome, bal.test17.f$outcome, bal.test18.f$outcome,
                                 bal.test19.f$outcome,  bal.test20.f$outcome,  
                                 bal.test23.f$outcome, bal.test24.f$outcome, bal.test25.f$outcome, bal.test26.f$outcome,
                                 bal.test27.f$outcome),
                       Estimate =   c(bal.test1.f$coefficients[3],bal.test2.f$coefficients[3],bal.test3.f$coefficients[3],bal.test4.f$coefficients[3],bal.test5.f$coefficients[3],
                                      bal.test6.f$coefficients[3],bal.test7.f$coefficients[3],bal.test8.f$coefficients[3],bal.test9.f$coefficients[3],bal.test10.f$coefficients[3],
                                      bal.test11.f$coefficients[3],bal.test12.f$coefficients[3],bal.test13.f$coefficients[3], bal.test14.f$coefficients[3],
                                      bal.test15.f$coefficients[3], bal.test16.f$coefficients[3], bal.test17.f$coefficients[3], bal.test18.f$coefficients[3],
                                      bal.test19.f$coefficients[3],  bal.test20.f$coefficients[3], 
                                      bal.test23.f$coefficients[3], bal.test24.f$coefficients[3], bal.test25.f$coefficients[3], bal.test26.f$coefficients[3],
                                      bal.test27.f$coefficients[3]),
                       pvalue =     c(bal.test1.f$p.value[3],bal.test2.f$p.value[3],bal.test3.f$p.value[3],bal.test4.f$p.value[3],bal.test5.f$p.value[3],
                                      bal.test6.f$p.value[3],bal.test7.f$p.value[3],bal.test8.f$p.value[3],bal.test9.f$p.value[3],bal.test10.f$p.value[3],
                                      bal.test11.f$p.value[3],bal.test12.f$p.value[3],bal.test13.f$p.value[3], bal.test14.f$p.value[3],
                                      bal.test15.f$p.value[3], bal.test16.f$p.value[3], bal.test17.f$p.value[3], bal.test18.f$p.value[3],
                                      bal.test19.f$p.value[3], bal.test20.f$p.value[3], 
                                      bal.test23.f$p.value[3], bal.test24.f$p.value[3], bal.test25.f$p.value[3], bal.test26.f$p.value[3],
                                      bal.test27.f$p.value[3]),
                       SE =         c(bal.test1.f$std.error[3],bal.test2.f$std.error[3],bal.test3.f$std.error[3],bal.test4.f$std.error[3],bal.test5.f$std.error[3],
                                      bal.test6.f$std.error[3],bal.test7.f$std.error[3],bal.test8.f$std.error[3],bal.test9.f$std.error[3],bal.test10.f$std.error[3],
                                      bal.test11.f$std.error[3],bal.test12.f$std.error[3],bal.test13.f$std.error[3], bal.test14.f$std.error[3],
                                      bal.test15.f$std.error[3], bal.test16.f$std.error[3], bal.test17.f$std.error[3], bal.test18.f$std.error[3],
                                      bal.test19.f$std.error[3], bal.test20.f$std.error[3], 
                                      bal.test23.f$std.error[3], bal.test24.f$std.error[3], bal.test25.f$std.error[3], bal.test26.f$std.error[3],
                                      bal.test27.f$std.error[3]))
balance1$Coefficient <- ifelse(balance1$Estimate>0, "Positive", "Negative")
balance1$Var <- factor(balance1$Var, ordered(balance1$Var))
balance1$Specification <- "Full"

balance.data <- rbind(balance,balance1)

# Replace variables
balance.data$Var=c('Landless Population','Banks','Unemployment','Land Gini',
                   'Income Gini','Rain Deviation (Annual)','Education HDI', 'Log (Security Budget)',
                   'Log (Social Spending)','Agricultural Income','NonAgricultural GDP','Extreme Poverty',
                   'Distance from Capital', 'Invasions', 'Reforms', 'Ratio farms <1ha/>100has', 'Percent rural', 'Municipal guard',
                   'Agricultural productivity', 'Cattle dependency',  'Coffee dependency',
                   'Left governor', 'Right governor', 'PSDB vote share', 'PT vote share',
                   'Landless Population','Banks','Unemployment','Land Gini',
                   'Income Gini','Rain Deviation (Annual)','Education HDI', 'Log (Security Budget)',
                   'Log (Social Spending)','Agricultural Income','NonAgricultural GDP','Extreme Poverty',
                   'Distance from Capital', 'Invasions', 'Reforms', 'Ratio farms <1ha/>100has', 'Percent rural', 'Municipal guard',
                   'Agricultural productivity', 'Cattle dependency', 'Coffee dependency',
                   'Left governor', 'Right governor', 'PSDB vote share', 'PT vote share')

# balance plot
ggplot(balance.data)+ 
  geom_point(aes(x = pvalue, y=Var, color=Specification, shape=Specification), 
             lwd = 2, position = position_dodge(width = 0.5))+
  geom_vline(xintercept = 0.05, colour = "red", lty = 2, lwd=0.3)+
  scale_color_manual(values=c("royalblue2", "tomato2"))+
  xlim(0,1)+
  ylab('')+xlab('p-value')+theme_bw()+
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(hjust = 0),
        axis.text.y = element_text(hjust = 0),
        legend.position="bottom")
 

ggsave(file = "appendix_figures/Figure_B1.pdf",
       width = 11, height = 12.5, dpi = 600)


#================================================
#                 Figure B.2
#               McCrary's Test
#================================================

pdf("appendix_figures/Figure_B2.pdf",
    width = 10, height = 7)
DCdensity(bra_all$prog, 0) 
abline(v=0)
dev.off()


#========================================================
#                 Table B.3
#   Effects of political parties on violence - covariates
#========================================================

# covariates
covspp1 <- cbind(bra_all$ruralpc, bra_all$distcap)
covspp2 <- cbind(bra_all$ruralpc, bra_all$homlag, bra_all$distcap)

# programmatic parties
summary(m1_1<-rdrobust(bra_all$hom_avg,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd", covs=NULL))

summary(m1_2<-rdrobust(bra_all$hom_avg,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd", covs=covspp1))

summary(m1_3<-rdrobust(bra_all$hom_avg,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd", covs=covspp2))

msummary(list(m1_1, m1_2, m1_3),
         stars = T,
         output = "latex") %>%
  add_header_above(c(" " = 1, "Homicide rates (Avg)" = 3)) %>%
  add_header_above(c("Panel A: Programmatic party)" = 4))

# non-programmatic parties
summary(npm1_1<-rdrobust(bra_all$hom_avg,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd", covs=NULL))

summary(npm1_2<-rdrobust(bra_all$hom_avg,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd", covs=covspp1))

summary(npm1_3<-rdrobust(bra_all$hom_avg,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd", covs=covspp2))

msummary(list(npm1_1, npm1_2, npm1_3),
         stars = T,
         output = "latex") %>%
  add_header_above(c(" " = 1, "Homicide rates (Avg)" = 3)) %>%
  add_header_above(c("Panel B: Non-Programmatic party)" = 4))

#================================================
#                 Table 
#       Disaggregating by political party
#================================================

## programmatic parties

# PT only
summary(pp_1<-rdrobust(bra_all$hom_avg,bra_all$prog_1,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd"))
# PT & PSDB 
summary(pp_2<-rdrobust(bra_all$hom_avg,bra_all$prog_2,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd"))
# PT, PSDB, PDT, & PP
summary(pp_3<-rdrobust(bra_all$hom_avg,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd"))

msummary(list(pp_1, pp_2, pp_3),
         stars = T,
         output = "latex") %>%
  add_header_above(c(" " = 1, "Homicide rates (Avg)" = 3)) %>%
  add_header_above(c("Panel A: Programmatic party)" = 4))

## non-programmatic parties

# PMDB only
summary(np_1<-rdrobust(bra_all$hom_avg,bra_all$nonprog_1,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd"))

# PTB, PPB, PL, PSB
summary(np_2<-rdrobust(bra_all$hom_avg,bra_all$nonprog_2,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd"))

# PMDB, PSB, PFL, PPB, PTB
summary(np_3<-rdrobust(bra_all$hom_avg,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd"))
#export
msummary(list(np_1, np_2, np_3),
         stars = T, output = "latex") %>%
  add_header_above(c(" " = 1, "Homicide rates (Avg)" = 3)) %>%
  add_header_above(c("Panel B: Non-Programmatic party)" = 4))



#================================================
#                    Table 
#           Effects of party over time
#================================================

## programmatic parties

summary(m2_1<-rdrobust(bra_all$hom0,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd"))
summary(m3_1<-rdrobust(bra_all$hom1,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd"))
summary(m4_1<-rdrobust(bra_all$hom2,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd"))
summary(m5_1<-rdrobust(bra_all$hom3,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd"))
summary(m6_1<-rdrobust(bra_all$hom4,bra_all$prog,c=0,p=1,  
                       all = T, vce = "hc0", kernel="tri",
                       bwselect = "mserd"))

#export
msummary(list(m2_1, m3_1, m4_1, m5_1, m6_1),
         stars = T, output = "latex") %>%
  add_header_above(c(" " = 1, "Homicide rates (Avg)" = 5)) %>%
  add_header_above(c("Panel A: Programmatic party)" = 6))

## non-programmatic parties
summary(npm2_1<-rdrobust(bra_all$hom0,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd"))
summary(npm3_1<-rdrobust(bra_all$hom1,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd"))
summary(npm4_1<-rdrobust(bra_all$hom2,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd"))
summary(npm5_1<-rdrobust(bra_all$hom3,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd"))
summary(npm6_1<-rdrobust(bra_all$hom4,bra_all$nonprog,c=0,p=1,  
                         all = T, vce = "hc0", kernel="tri",
                         bwselect = "mserd"))

#export
msummary(list(npm2_1, npm3_1, npm4_1, npm5_1, npm6_1),
         stars = T, output = "latex") %>%
  add_header_above(c(" " = 1, "Homicide rates (Avg)" = 5)) %>%
  add_header_above(c("Panel B: Non-Programmatic party)" = 6))



#================================================
#                    Figure 
#       Placebo test - trying different cutoffs
#================================================

# create placebo cutoffs
cutoff <- sort(c(seq(-.35,.35,by=.02),0))

# matrix to store homicide estimates
store_p_hom <- as.data.frame(matrix(NA,37,3))
colnames(store_p_hom)<-c('coef','se','pv')


# Initiate loop for homicides #
for (i in 1:length(cutoff)) {
  # Estimate model
  summary(hom_loop <- rdrobust(y=bra_all$hom_avg, x = bra_all$prog, c = cutoff[i], p=1, vce = "hc0", 
                               bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_p_hom[[i,1]] <- hom_loop$coef[3]
  store_p_hom[[i,2]] <- hom_loop$se[3]
  store_p_hom[[i,3]] <- hom_loop$pv[3]
  store_p_hom$c <- cutoff
  store_p_hom$c<-round(store_p_hom$c, digits = 2)
  
}



store_p_hom %>%
  ggplot()+ 
  geom_errorbar(aes(x = factor(c), ymin = coef - se*1.96,
                    ymax = coef + se*1.96), color="red", lwd = .93, width=.3)+
  geom_point(aes(x = factor(c), y = coef), lwd = 3, shape = 20)+
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)+
  ylim(-15,15) + xlab('') + ylab('') + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45, hjust = 1, size=15))


ggsave(file = "appendix_figures/Figure_B4.pdf",
       width = 10, height = 8, dpi = 600)



#================================================
#                    Figure 
#       RD estimates using different bandwidths
#================================================


# create bandwidths
bandw <- seq(.05,.21,by=.005)

# matrix to store homicide estimates
store_bw1 <- as.data.frame(matrix(NA,33,3)); store_bw1.np<- as.data.frame(matrix(NA,33,3))
store_bw2 <- as.data.frame(matrix(NA,33,3)); store_bw2.np<- as.data.frame(matrix(NA,33,3))
store_bw3 <- as.data.frame(matrix(NA,33,3)); store_bw3.np<- as.data.frame(matrix(NA,33,3))
store_bw4 <- as.data.frame(matrix(NA,33,3)); store_bw4.np<- as.data.frame(matrix(NA,33,3))
store_bw5 <- as.data.frame(matrix(NA,33,3)); store_bw5.np<- as.data.frame(matrix(NA,33,3))
store_bw6 <- as.data.frame(matrix(NA,33,3)); store_bw6.np<- as.data.frame(matrix(NA,33,3))


# Initiate loop for homicides t
for (i in 1:length(bandw)) {
  # Estimate model
  summary(bw_loop1 <- rdrobust(y=bra_all$hom0, x = bra_all$prog, c =0, p=1,h=bandw[i],vce = "hc0", 
                               bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw1[[i,1]] <- bw_loop1$coef[3]
  store_bw1[[i,2]] <- bw_loop1$se[3]
  store_bw1[[i,3]] <- bw_loop1$pv[3]
  store_bw1$c <- bandw
  store_bw1$Party <- "Programmatic"
  store_bw1$c<-round(store_bw1$c, digits = 4)
  colnames(store_bw1)<-c('coef','se','pv', 'c', 'Party')
  
  # Estimate model
  summary(bw_loop1.np <- rdrobust(y=bra_all$hom0, x = bra_all$nonprog, c =0, p=1,h=bandw[i],vce = "hc0", 
                                  bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw1.np[[i,1]] <- bw_loop1.np$coef[3]
  store_bw1.np[[i,2]] <- bw_loop1.np$se[3]
  store_bw1.np[[i,3]] <- bw_loop1.np$pv[3]
  store_bw1.np$c <- bandw
  store_bw1.np$Party <- "Non-programmatic"
  store_bw1.np$c<-round(store_bw1.np$c, digits = 4)
  colnames(store_bw1.np)<-c('coef','se','pv', 'c', 'Party')
  store1<-rbind(store_bw1,store_bw1.np)  
}


# Initiate loop for homicides t+1
for (i in 1:length(bandw)) {
  # Estimate model
  summary(bw_loop2 <- rdrobust(y=bra_all$hom1, x = bra_all$prog, c =0, p=1,h=bandw[i],vce = "hc0", 
                               bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw2[[i,1]] <- bw_loop2$coef[3]
  store_bw2[[i,2]] <- bw_loop2$se[3]
  store_bw2[[i,3]] <- bw_loop2$pv[3]
  store_bw2$c <- bandw
  store_bw2$Party <- "Programmatic"
  store_bw2$c<-round(store_bw2$c, digits = 4)
  colnames(store_bw2)<-c('coef','se','pv', 'c', 'Party')
  
  # Estimate model
  summary(bw_loop2.np <- rdrobust(y=bra_all$hom1, x = bra_all$nonprog, c =0, p=1,h=bandw[i],vce = "hc0", 
                                  bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw2.np[[i,1]] <- bw_loop2.np$coef[3]
  store_bw2.np[[i,2]] <- bw_loop2.np$se[3]
  store_bw2.np[[i,3]] <- bw_loop2.np$pv[3]
  store_bw2.np$c <- bandw
  store_bw2.np$Party <- "Non-programmatic"
  store_bw2.np$c<-round(store_bw2.np$c, digits = 4)
  colnames(store_bw2.np)<-c('coef','se','pv', 'c', 'Party')
  store2<-rbind(store_bw2,store_bw2.np)  
}



# Initiate loop for homicides t+2
for (i in 1:length(bandw)) {
  # Estimate model
  summary(bw_loop3 <- rdrobust(y=bra_all$hom2, x = bra_all$prog, c =0, p=1,h=bandw[i],vce = "hc0", 
                               bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw3[[i,1]] <- bw_loop3$coef[3]
  store_bw3[[i,2]] <- bw_loop3$se[3]
  store_bw3[[i,3]] <- bw_loop3$pv[3]
  store_bw3$c <- bandw
  store_bw3$Party <- "Programmatic"
  store_bw3$c<-round(store_bw3$c, digits = 4)
  colnames(store_bw3)<-c('coef','se','pv', 'c', 'Party')
  
  # Estimate model
  summary(bw_loop3.np <- rdrobust(y=bra_all$hom2, x = bra_all$nonprog, c =0, p=1,h=bandw[i],vce = "hc0", 
                                  bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw3.np[[i,1]] <- bw_loop3.np$coef[3]
  store_bw3.np[[i,2]] <- bw_loop3.np$se[3]
  store_bw3.np[[i,3]] <- bw_loop3.np$pv[3]
  store_bw3.np$c <- bandw
  store_bw3.np$Party <- "Non-programmatic"
  store_bw3.np$c<-round(store_bw3.np$c, digits = 4)
  colnames(store_bw3.np)<-c('coef','se','pv', 'c', 'Party')
  store3<-rbind(store_bw3,store_bw3.np)  
}

# Initiate loop for homicides t+3
for (i in 1:length(bandw)) {
  # Estimate model
  summary(bw_loop4 <- rdrobust(y=bra_all$hom3, x = bra_all$prog, c =0, p=1,h=bandw[i],vce = "hc0", 
                               bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw4[[i,1]] <- bw_loop4$coef[3]
  store_bw4[[i,2]] <- bw_loop4$se[3]
  store_bw4[[i,3]] <- bw_loop4$pv[3]
  store_bw4$c <- bandw
  store_bw4$Party <- "Programmatic"
  store_bw4$c<-round(store_bw4$c, digits = 4)
  colnames(store_bw4)<-c('coef','se','pv', 'c', 'Party')
  
  # Estimate model
  summary(bw_loop4.np <- rdrobust(y=bra_all$hom3, x = bra_all$nonprog, c =0, p=1,h=bandw[i],vce = "hc0", 
                                  bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw4.np[[i,1]] <- bw_loop4.np$coef[3]
  store_bw4.np[[i,2]] <- bw_loop4.np$se[3]
  store_bw4.np[[i,3]] <- bw_loop4.np$pv[3]
  store_bw4.np$c <- bandw
  store_bw4.np$Party <- "Non-programmatic"
  store_bw4.np$c<-round(store_bw4.np$c, digits = 4)
  colnames(store_bw4.np)<-c('coef','se','pv', 'c', 'Party')
  store4<-rbind(store_bw4,store_bw4.np)  
}



# Initiate loop for homicides t+4
for (i in 1:length(bandw)) {
  # Estimate model
  summary(bw_loop5 <- rdrobust(y=bra_all$hom4, x = bra_all$prog, c =0, p=1,h=bandw[i],vce = "hc0", 
                               bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw5[[i,1]] <- bw_loop5$coef[3]
  store_bw5[[i,2]] <- bw_loop5$se[3]
  store_bw5[[i,3]] <- bw_loop5$pv[3]
  store_bw5$c <- bandw
  store_bw5$Party <- "Programmatic"
  store_bw5$c<-round(store_bw5$c, digits = 4)
  colnames(store_bw5)<-c('coef','se','pv', 'c', 'Party')
  
  # Estimate model
  summary(bw_loop5.np <- rdrobust(y=bra_all$hom4, x = bra_all$nonprog, c =0, p=1,h=bandw[i],vce = "hc0", 
                                  bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw5.np[[i,1]] <- bw_loop5.np$coef[3]
  store_bw5.np[[i,2]] <- bw_loop5.np$se[3]
  store_bw5.np[[i,3]] <- bw_loop5.np$pv[3]
  store_bw5.np$c <- bandw
  store_bw5.np$Party <- "Non-programmatic"
  store_bw5.np$c<-round(store_bw5.np$c, digits = 4)
  colnames(store_bw5.np)<-c('coef','se','pv', 'c', 'Party')
  store5<-rbind(store_bw5,store_bw5.np)  
}



# Initiate loop for avg homicides  
for (i in 1:length(bandw)) {
  # Estimate model
  summary(bw_loop6 <- rdrobust(y=bra_all$hom_avg, x = bra_all$prog, c =0, p=1,h=bandw[i],vce = "hc0", 
                               bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw6[[i,1]] <- bw_loop6$coef[3]
  store_bw6[[i,2]] <- bw_loop6$se[3]
  store_bw6[[i,3]] <- bw_loop6$pv[3]
  store_bw6$c <- bandw
  store_bw6$Party <- "Programmatic"
  store_bw6$c<-round(store_bw6$c, digits = 4)
  colnames(store_bw6)<-c('coef','se','pv', 'c', 'Party')
  
  # Estimate model
  summary(bw_loop6.np <- rdrobust(y=bra_all$hom_avg, x = bra_all$nonprog, c =0, p=1,h=bandw[i],vce = "hc0", 
                                  bwselect = "mserd", kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw6.np[[i,1]] <- bw_loop6.np$coef[3]
  store_bw6.np[[i,2]] <- bw_loop6.np$se[3]
  store_bw6.np[[i,3]] <- bw_loop6.np$pv[3]
  store_bw6.np$c <- bandw
  store_bw6.np$Party <- "Non-programmatic"
  store_bw6.np$c<-round(store_bw6.np$c, digits = 4)
  colnames(store_bw6.np)<-c('coef','se','pv', 'c', 'Party')
  store6<-rbind(store_bw6,store_bw6.np)  
}



## plots 
# homicide rates - t
store1 %>%
  ggplot()+ 
  geom_errorbar(aes(x = factor(c), ymin = coef - se*1.96,
                    ymax = coef + se*1.96, color=Party), lwd = 1, width=.3,
                position=position_dodge(width=0.9))+
  geom_point(aes(x = factor(c), y = coef, color=Party), lwd = 2, shape = 20,
             position=position_dodge(width=0.9))+
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)+
  ylim(-15,25)+ xlab('') + ylab('') + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 11),
        legend.position="bottom")+
  scale_color_manual(values=c("#999999", "green4"))

ggsave(file = "appendix_figures/Figure_B5a.pdf",
       width = 8, height = 7, dpi = 600)

# homicide rates - t+1
store2 %>%
  ggplot()+ 
  geom_errorbar(aes(x = factor(c), ymin = coef - se*1.96,
                    ymax = coef + se*1.96, color=Party), lwd =1, width=.3,
                position=position_dodge(width=0.9))+
  geom_point(aes(x = factor(c), y = coef, color=Party), lwd = 2, shape = 20,
             position=position_dodge(width=0.9))+
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)+
  ylim(-15,25)+ xlab('') + ylab('') + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 11),
        legend.position="bottom")+
  scale_color_manual(values=c("#999999", "green4"))


ggsave(file = "appendix_figures/Figure_B5b.pdf",
       width = 8, height = 7, dpi = 600)

# homicide rates - t+2
store3 %>%
  ggplot()+ 
  geom_errorbar(aes(x = factor(c), ymin = coef - se*1.96,
                    ymax = coef + se*1.96, color=Party), lwd = 1, width=.3,
                position=position_dodge(width=0.9))+
  geom_point(aes(x = factor(c), y = coef, color=Party), lwd = 2, shape = 20,
             position=position_dodge(width=0.9))+
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)+
  ylim(-15,25)+ xlab('') + ylab('') + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 11),
        legend.position="bottom")+
  scale_color_manual(values=c("#999999", "green4"))


ggsave(file = "appendix_figures/Figure_B5c.pdf",
       width = 8, height = 7, dpi = 600)


# homicide rates - t+3
store4 %>%
  ggplot()+ 
  geom_errorbar(aes(x = factor(c), ymin = coef - se*1.96,
                    ymax = coef + se*1.96, color=Party), lwd = 1, width=.3,
                position=position_dodge(width=0.9))+
  geom_point(aes(x = factor(c), y = coef, color=Party), lwd = 2, shape = 20,
             position=position_dodge(width=0.9))+
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)+
  ylim(-15,25)+ xlab('') + ylab('') + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 11),
        legend.position="bottom")+
  scale_color_manual(values=c("#999999", "green4"))

ggsave(file = "appendix_figures/Figure_B5d.pdf",
       width = 8, height = 7, dpi = 600)

# homicide rates - t+4
store5 %>%
  ggplot()+ 
  geom_errorbar(aes(x = factor(c), ymin = coef - se*1.96,
                    ymax = coef + se*1.96, color=Party), lwd =1, width=.3,
                position=position_dodge(width=0.9))+
  geom_point(aes(x = factor(c), y = coef, color=Party), lwd = 2, shape = 20,
             position=position_dodge(width=0.9))+
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)+
  ylim(-15,25)+ xlab('') + ylab('') + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 11),
        legend.position="bottom")+
  scale_color_manual(values=c("#999999", "green4"))

ggsave(file = "appendix_figures/Figure_B5e.pdf",
       width = 8, height = 7, dpi = 600)

# homicide rates - entire mayoral term
store6 %>%
  ggplot()+ 
  geom_errorbar(aes(x = factor(c), ymin = coef - se*1.96,
                    ymax = coef + se*1.96, color=Party), lwd = 1, width=.3,
                position=position_dodge(width=0.9))+
  geom_point(aes(x = factor(c), y = coef, color=Party), lwd = 2, shape = 20,
             position=position_dodge(width=0.9))+
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)+
  ylim(-15,25)+ xlab('') + ylab('') + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 11),
        legend.position="bottom")+
  scale_color_manual(values=c("#999999", "green4"))

ggsave(file = "appendix_figures/Figure_B5f.pdf",
       width = 8, height = 7, dpi = 600)

